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ABSTRACT 

This paper presents a continuation of our efforts to numerically study accretion disks that are 
misaligned (tilted) with respect to the rotation axis of a Kerr black hole. Here we present results 
of a global numerical simulation which fully incorporates the effects of the black hole spacetime as 
well as magnetorotational turbulence that is the primary source of angular momentum transport 
in the flow. This simulation shows dramatic differences from comparable simulations of untiltcd 
disks. Accretion onto the hole occurs predominantly through two opposing plunging streams that 
start from high latitudes with respect to both the black-hole and disk midplanes. This is due 
to the aspherical nature of the gravitational spacetime around the rotating black hole. These 
plunging streams start from a larger radius than would be expected for an untilted disk. In this 
regard the tilted black hole effectively acts like an untilted black hole of lesser spin. Throughout 
the duration of the simulation, the main body of the disk remains tilted with respect to the 
symmetry plane of the black hole; thus there is no indication of a Bardeen-Petterson effect in 
the disk at large. The torque of the black hole instead principally causes a global precession of 
the main disk body. In this simulation the precession has a frequency of 3(M Q /M) Hz, a value 
consistent with many observed low-frequency quasi-periodic oscillations. However, this value is 
strongly dependent on the size of the disk, so this frequency may be expected to vary over a large 
range. 

Subject headings: accretion, accretion disks — black hole physics — galaxies: active — MHD 
relativity — X-rays: stars 



Introduction 



Black-hole accretion has long been postulated to power the energetic emissions seen from quasars, 
active galactic nuclei (AGN), and many galactic X-ray sources: there is n ow ample observational evidence 



to support such claims (e.g. iKrolild Il999t iMcClintock fc Remillardl I2005T ) . Black- hole accretion flows are 



also of interest as laboratories to test predictions of general relativity. However, the nature of such flows 
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is complex, involving time-dependent, multi-dimensional dynamics with generically little symmetry. Hence 
numerical simulations play an integral role in advancing our understanding. 

Many simulations o f black- hole accretion flows have been carried ou t over the past three decades, both in 
the hydrodynamic (e.g. IWilsonlll97 2; Hawlc v et al.lll 984: lHawlev|Jl99ll) and magnetohydrodynamic (MHD) 
(e.g. iKoide et al.lll999HGammie et al.ll2003l 



De Villiers fc HawlevlEooibT) regimes. A common assumption in 



nearly all of the work to date has been that the symmetry plane of the central black hole is aligned with the 
midplane of the accretion flow, at least in some averaged sense. However, there is compelling ob servational 
evidence in several black-hole X-ray binaries (BHBs), e.g. GRO J1655-40 (jOrosz fe Bailvn 1997) and XTE 



J1550-564 (Hannikainenetal. 2001: Orosz et al 



NGC 1068 (|Caproni et alJl2006h . and NGC 4258 ( Caproni et al.ll2007l) . suggesting that misaligned (or tilted) 



black holes may be common (see also iMaccarondbOO^ This claim relies on the observation of relativistic 
bipolar jets (thought to be aligned with the spin axis of the black hole) that are not perpendicular to the 
plane of the accretion disk observed at large scales. 

There are also compelling theoretical arguments that many black holes should be tilted. First, the 
formation avenues for many black-hole - disk systems favor, or at least allow for, a tilted configuration 
([Fragile et al.ll200ll ). In stellar mass binaries, the orientation of the outer disk is fixed by the binary orbit, 
whereas the orientation of the black hole is determined by how it became part of the system, whether through 
a supernova explosion or multi-body interaction. If the black hole formed from a member of a preexisting 
binary through a supernova, then the black hole could be tilted if the explosion were asymmetric. If the 
black hole joined the binary through multi-body interactions, such as binary capture or replacement, then 
there would have been no preexisting symmetry, so the resulting system would nearly always harbor a tilted 
black hole. This same argument can be extended to AGN in which merger events reorient the central black 
hole or its fuel supply and result in repeated tilted configurations. 

If an accretion disk is misaligned or tilted, it will be subject to Lense-Thirring precession. For an ideal 
test particle in a slightly tilted orbit at a radius r around a black hole of mass M and specific angular 
momentum a, this precession occurs at an angular frequency Q^t ~ 2aM/r 3 . Close to the black hole, this is 
comparable to the orbital angular frequency £1 = (A//r 3 ) 1 / 2 /[l + a(M/r 3 ) 1 / 2 ] ss ^Kep- However, because of 
its strong radial dependence, Lense-Thirring precession becomes much weaker far from the hole. Therefore, 
a disk will experience a differential precession that will tend to twist and warp it. 

A warping disturbance can be communicated through a disk in either a diffusive or wave- like manner. 
In the diffusive case, the warping is limited by secular (i.e. "viscous") responses within the disk. In such a 
case, Lense-Thirring precess i on is expected to domin ate out to a unique, nearly constant transition radius 



20021). and AGN. e. g. NGC 3079 dKondratko et~al1l2005h . 



( Bardeen fc Petterson 1975 : Kumar fc Pringlel 1985 ). inside of which the disk is expected to be flat and 



aligned with the black-hole midplane, and outside of which the disk is also expected to be flat but in a 
plane determined by the angular momentum vector of the gas reservoir. This is what we term a "Bardeen- 
Petterson" configuration. Interestingly, data for the two black-hole X-ray binaries previously mentioned 
are best fit by disk components with inclinations that differ from their binary measure ments. The best-fi t 
inclinations are more consistent with inclination const raints derived from t he radio jets (jDavis et al.l l2006). 
possibly suggesting Bardeen-Petterson configurations. I Caproni et al. ( 20061 ) also claim that the observations 
of NGC 1068 are consistent with the predictions of the Bardeen-Pett erson effect. Confir mation could come 
through observations of relativistically broadened reflection features ([Fragile et al.ll2005T ). 



The Bar deen-Petterson result is ex pected to apply for Keplerian disks whenever the dimensionless stress 
parameter a ( Shakura fc Sunvaevlll973 ) is larger than the ratio of the disk semi-thickness H to the radius r at 
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all radii. Given that a is usually considered to be significantly less than one, this implies very geometrically 
"thin" disks. Unfortunately, current computational limitations prevent us from conducting global simulations 
of disks that are this thin. On the other hand, the Bardeen-Petterson regime may not be that common in 
real disks. Neglecting relativistic correction factors, the innermost, radiation pressure and electron scattering 
dominated portions of radiatively efficient accretion disks satisfy 

r \ L Edd/ \r G J 

where e ~ 0.1 is the radiative efficiency, i/LEdd is the luminosity in units of Eddington, and tq = GM/c 2 
is the gravitational radius. Note that equation (JTJ) is independent of whether the stress is chosen to be 
proportional to gas pressure, radiation pressure, or some combination of the two. We therefore conclude 
that the Bardeen-Petterson regime will be relevant in radiatively efficient disks near the black hole only for 
very small Eddington ratios L/L^dd S a£ << !• Moreover, radiatively less efficient, geometrically slim and 
thick flows will clearly not be in the Bardeen-Petterson regime. 

Global simulations of tilted disks that have H/r > a are computationally feasible. In this regime Lense- 



Thirring p recession is expecte d to pr oduce warps that propagate in a wave- like manner (|Papaloizou fe Lin 



1995 ). In Fragile fc Anninosl ( 2005 ) we presented results from the first fully general relativistic three- 



dimensional hydrodynamic numerical studies of tilted thick-disk accretion onto rapidly rotating (Kerr) black 
holes. We found that, although Lense-Thirring precession did cause the disk to warp, the warping only 
occurred inside a radius in the disk at which the precession time became comparable to other dynamical 
timescales, primarily the azimuthal sound-crossing time. After the differential warping ended and the evo- 
lution became quasi-static, the disks underwent near solid-body precession at rates consistent with some 
low- frequency quasi-periodic oscillations (QPOs). 

In this paper we extend the results of Fragile fc Anninosl ( 2005 ) to include magnetic fields. The 



in- 



clusion of magnetic fields is important because it is now widely believed that local stresses within black- 
hole accretion disks are generated by turbulence that results from the magnetorotational instability (MRI; 



Balbus fc Hawlevlll99lf ). Here we report on our first global general relativistic MHD (GRMHD) simulation 
of a tilted accretion disk around a moderately rapidly rotating black hole (a/M = 0.9). The simulation is 
initialized starting from the analytic solution for an axisymmetric torus around a rotating black hole. A 
weak poloidal magnetic field is added to the torus to seed the MRI. After the torus is initialized, the black 
hole is tilted by an angle /3q = 15° relative to the disk through a transformation of the metric. The system is 
then allowed to evolve. This paper reports the results as follows: In $2]we describe the numerical procedures 
used in this GRMHD simulation. In <j3]we present the results of this simulation. In <|4]we summarize our 
findings and draw conclusions. 



Numerical Methods 



This work is carried out using the Cosmos++ as trophysical magnetohydro dynamics code (jAnninos et al 



20051 ). Similar to our predecessor code Cosmos (jAnninos fc Fragile! 120031 ). Cosmos++ includes several 



schemes for solving the GRMHD equations. The fluid equations can be solved using a traditional artifi- 
cial viscosity scheme, non-oscillatory central difference methods, or a new hybrid dual energy (internal and 
total) method. For this work, we use the artificial viscosity formulation, mainly because of its speed and 
robustness. With the magnetic fields we solve the induction equation in an advection-split form and apply a 
hyperbolic divergence cleanser to maintain an approximately divergence-free magnetic field. For clarity and 
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notation sake, we present the full evolution equations for mass, internal energy, momentum, and magnetic 
induction as solved in this work. Throughout this paper we use units where G = c = 1 and the metric sig- 
nature is (—,+,+,+). We use the standard notation in which four- and three-dimensional tensor quantities 
are represented by Greek and Latin indices, respectively. 

The evolution equations are 

d t D + d t (DV l ) = 0, (2) 
dtE + d^EV 1 ) = -Pd t W-(P + Q)d l (WV i ) , (3) 

+ (f^ - ^ B " BU ) d ^ V=S dj (P + Pb + Q) , (4) 
d t B j + d i (B j V i ) = tf'),\' J ■ g ij <),r . (5) 



.2 



d^ + c^B 1 = — (6) 
c p 

where is the 4-metric, g is the 4-metric determinant, W = ^—gvP is the relativistic boost factor, 
D = W p is the generalized fluid density, V 1 — u l /u° is the transport velocity, it M = g^ LV u u is the fluid 
4- velocity, = W(ph + 2Pb)u i1 is the covariant momentum density, E = We — Wpe is the generalized 
internal energy density, P is the fluid pressure, Q is the artificial viscosity used for shock capturing, and Ch 
and c p are coefficients to determine the strength of the hyperbolic and parabolic pieces of the divergence 
cleanser. There are two representations of the magnetic field in these equations: is the rest frame 
magnetic induction used in defining the stress tensor 

T" v = (ph + 2P B ) u^u v + (P + P B ) - —B tl B lJ (7) 

47T 

and 

B" = W{B» - B°V^) (8) 

is the divergence- free (dB 1 jdx % = 0), spatial (B° = 0) representation of the field. The time component of 
the magnetic field B° is recovered from the orthogonality condition B^u^ = 

B° = -j (goiB* + g^V*) . (9) 

The relativistic enthalpy is 

TP O 

k=1+ (T-Wp + P> (10) 
where we have assumed an equation of state of the form P = (T — l)pe. Finally, Ph = ||i?|| 2 /87r = 
g f _ l uB f *B l/ /8tt is the magnetic pressure. We use the scalar Q from lAnninos et aL ( 20051 ) with k q — 2.0 and 



ki = 0.3. We fix the divergence cleanser coefficients to be = c c flAa; m i n /At and c 2 = c/,., where c c fi = 0.7 
is the Courant coefficient, Ax m i n is the minimum covariant zone length, and At is the evolution timestep. 
For simplicity, we hold the timestep fixed at At = c c flAx m i n throughout the simulation. 

These GRMHD equations are evolved in a "tilted" Kerr-Schild polar coordinate system (t, r, d, cp). This 
coordinate system is related to the usual (untilted) Kerr-Schild coordinates (t, r, 9, </>) through a simple 
rotation about the y-axis by an angle Po, such that 

sin d cos (p \ I cos /3o — sin A) \ / sin 9 cos < 

sintfsin^ = 1 sin^sin^ | . (11) 

cos-d I V sin/3n COS/3n / V cos( 
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The full tilted metric terms are provided in Fragile fc Anninosl ( 2005 ) [see also Fragile fc Anninosl ( 2007 )]. 



The comput ational advantages of t he " horizon- adapted" Kerr-S child form of the Kerr metric were first 
described in IPapadopoulos &: Font ( 19981 ) and Font et al. ( 19981 ). The primary advantage is that, unlike 



Boyer-Lindquist coordinates, there are no singularities in the metric terms at the event horizon, so the 
computational mesh can extend into the hole's interior. In principle, this should keep the inner boundary 
causally disconnected from the flow, although numerically there is still some communication. 

The simulation is carried out on a spherical polar mesh with nested resolution layers. The base grid 
contains 32 3 mesh zones and covers the full 4ir steradians. Varying levels of refinement are added on top of this 
base layer; each refinement level doubles the resolution relative to the previous layer. The main simulation, 
referenced as Model 915h, has two levels of refinement, thus achieving a peak resolution equivalent to a 
128 3 simulation. For comparison we also discuss results from an equivalent untilted simulation (Model 90h) 
with the same resolution. As an argument that our results are reasonably well converged, we also include 
results from two other tilted simulations: one with a single refinement layer and an equivalent resolution 
of 64 3 (Model 915m) and another that starts from a base grid of 24 x 24 x 32 and adds three layers of 
refinement for an equivalent resolution of 192 x 192 x 256 (Model 915vh). The evolution times for these 
simulations differ as described below. In all cases, the full refinement covers the region r m j n < r < r" max , 
0.075tt = <8 X < <Q < d 2 = 0.925tt, < ip < 2tt, where r min = 0.98r BH = 1.41r G and r max = 120r G are 
the inner and outer boundaries of the grid, respectively, and tbh = 1.43r G is the black-hole horizon radius. 
The primary motivation for using a nested grid is to allow us to maintain a reasonable Courant-limited 
timestep without sacrificing any spatial resolution within the disk nor completely excluding the region near 
the pole. The gain in computational efficiency is significant since, for a polar mesh, the timestep scales as 
At r m ; n sin$ m i n A</?. By underresolving the polar region, we gain by increasing both i9 m ; n and Aip. With 
2 levels of refinement, we are able to use a timestep that is a factor of 11.8 larger than what we could use if 
our most refined layer extended all the way to the pole. The main drawback of this approach is that we are 
unable to resolve the region in which jets are expected to form. 

In the radial direction we use a logarithmic coordinate of the form r\ = 1.0 + ln(r/rBn)- The spatial 
resolution near the black-hole horizon is Ar « 0.05r G ; near the initial pressure maximum of the torus, the 
resolution is Ar w O.brc- Both are considerably smaller than the initial characteristic MRI wavelength 
Amri = 2ttva/Q ~ 2.5r G . This also gives us a large number of zones inside the plunging region. In the 
angular direction, in addition to the nested grids, we use a concentrated latitude coordinate xi of the form 
■d = X2 + |(1 — h) sin(2a;2) with h = 0.5, which concentrates resolution toward the midplane of the disk. As 
a result r cen terA?9 = 0.3r G near the midplane while it is a factor of ~ 3 larger for the fully refined zones near 
the pole. The grid used in Models 915h and 90h is shown in Figure [1] 

Since we cover the full 4tt steradians, the only "external" boundaries are the inner and outer radial 
boundaries, where we apply outflow conditions: Fluid variables are set the same in the external boundary 
zone as in the neighboring internal zone, except for velocity, which is chosen to satisfy 

yr _ ( ^int when V r points off the grid , , . 

cxt \ — V£ t when V r points onto the grid . 

In the azimuthal direction we apply periodic boundaries at ip = and 2ir. Since CosmosH — h is a zone-centered 
code, we do not have to treat the pole ($ = or tt) directly. Instead unboosted scalar quantities, such as the 
gas pressure P, in the "ghost" zones across the pole are filled with real data from the corresponding zone 
located 180° away in azimuth. Unboosted vector quantities, such as velocity V 1 , are similarly filled with 
data from appropriate real zones, albeit with the signs reversed for the $ and ip components to maintain a 
consistent sense of coordinate direction across the pole. Boosted quantities, since they contain the metric 
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Fig. 1. — Plot of the grid geometry used for the main simulation (Model 915h) 
in the symmetry plane of the grid, while the black hole is not. 



The initial torus is aligned 
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determinant y/—g, are reflected across the pole so they extr apolate to zero there. This treatment differs 
from the pure reflecting boundaries used in other works (e.g. De Villiers et al. 2003 : McKinnev 20061 ) in its 
treatment of the unboosted variables. For untiltcd black holes the difference is relatively minor. However, 
for tilted black holes, our approach makes the pole more transparent to the fluid. 

We initialize th ese simulations sta rting from the analytic solution for an axisymmetric torus around a 
rotating black hole ( Chakrabarti 19851) . To provide a link with an untilted model already in the literature, 
we start with identical torus conditions as model KDP of foe Villiers et al l |2003l ). which is the relativistic 
analog of model GT4 of lHawlevi (|2000f h In our initialization, the torus is defined by: the black- hole spacetime, 
specifically the spin of the black hole; the inner radius of the torus the radius of the pressure maximum of 
the torus r cen ter; and the power-law exponent q used in defining the specific angular momentum distribution, 



-u$/u t = kA 



2-9 



(13) 



As in model KDP, a/M = 0.9, r in = !5rQ, r ccntcr = 25r<3, and q = 1.68. Knowledge of r ccntor leads directly 
to a determination of £ cen ter by setting it equal to the geodesic value at that radius. The numerical value of 
k comes directly from the choice of q and the determination of A ccntcr , where 



A 2 ~ lg^ + £ 2 g t <t, 



(14) 



Finally, having chosen we can obtain itj n = Ut(ri n ), the surface binding energy of the torus, from 



9" 



2£g tct > 



"9" 



2003) 



The solution of the torus variables can now be specified. The internal energy of the torus is (|De Villiers et al 

V*inj \£in) 



<r,6) = - 



where L. 



u t {r,0)f{l{r,6))_ 
is the specific angular momentum of the fluid at the surface and 

f(£)= l-fc 2 /"f 



(15) 



(16) 



where n — 2 — q and a = (2n — 2)/n. Assuming an isentropic equation of state P — pe(T — 1) = np r , the 
density is given by p — [e(T — l)/^] 1 ^ 1 " 1 . As in model KDP, we take T = 5/3 and n = 0.01 (arbitrary 
units). Finally, the angular velocity of the fluid is specified by 



9tct> + ■ 



9 14> 



(17) 



The dependence of A on I in equation ([14]) for Kerr black holes means that the solution requires an 
iterative procedure. However, we can get an approximate solution by taking the Schwarzschild form (i.e. 
ignoring g t<j> ) 

p2 _ _gw_ _ (18) 
gtt 

The error introduced by doing so is small and only affects the initial torus configuration, which will already 
be unstable to the MRI due to the seed magnetic fields being added. Thus, this slightly si mplified treatment 



has no real consequence for the evolution. We note that the same procedure is followed in lDe Villiers et al 
d2003h . 
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Once the torus is constructed, it is seeded with a weak magnetic field in the form of poloidal loops along 
the iso baric contours within the torus. The initial magnetic field vector potential is (|De Villiers fc Hawlev 
2003ah 

KP ~ Pent) for p > p cut , 
for p < p cut . 



— 



(19) 



The non-zero spatial magnetic field components are then B r = —d^A v and = d r A v . The parameter 
Pcut = 0.5 * p max ,o is used to keep the field a suitable distance inside the surface of the torus, where p max ,o is 
the initial density maximum within the torus. Using the constant b in equation (|19p . the field is normalized 
such that initially @ m!l v = PJ Pr > /3mag,o = 10 throughout the torus. This initialization is slightly different 
than lDe Villiers fe Hawlevl (|2003bl ). who use a volume integrated f3 mag to set the field strength; the difference 
is such that /3 mag .o = 100 in their work is roughly comparable to /3 m ag,o = 10 here. 

In the background regio n not specified by t he torus solution, we set up a rarefied non-magnetic plasma 
accreting into the black hole (|Komissarovll2006l ) . The density and pressure have the form 



10 



-3 



Pmax,0 exp 



(— ) 

\ ^center / 



P = np v 



The radial velocity has the form 



yr 



fJ 1 



1 



(20) 



(21) 



This introduces inflow through the horizon without creating large veloci ty jumps at the t o rus su rface. This 
background is initially more dense than the static background used by iDe Villiers et al.l (|20031 ). However, 
since this background reservoir is not replenished at the outer boundary, it is rapidly depleted and has virtu- 
ally no long-term dynamical impact on the problem. Numerical floors are placed on p and e at approximately 
1CP 10 and 10~ 16 of their initial maxima, respectively. These floors are very seldom applied once the initial 
background is replaced by evolved disk material. 

The final step of the initialization is to tilt the black hole by an angle /3o = 15 ° relative to the disk (and 
the grid) by transforming the Kerr metric. The full transformation is provided in Fragile fc Anninos ( 2005h 
[see also iFragile fc Anninod (|2007T )]. Thus, while the torus is responding to the action of the MRI, it will 
also experience a gravitomagnetic torque from the tilted black hole. 



3. Results 

In the main simulation (915h) the torus is evolved for a total of 10 orbital periods (10i O rt>) as measured 
at r = r cen ter) which corresponds to ~ 350 orbits near nsco — 2.32r , <3, the coordinate radius of the innermost 
stable circular orbit (for prograde orbits in the symmetry plane of the black hole) . The very high resolution 
simulation (915vh) is only run for half as long (5£ rb), while the lower resolution simulation (915m) is run 
for twice as long (20i or b). Figure [2] shows snapshots of the disk from Model 915h at times t = 0, 1, 2, 4, 7, 
and 10t or b- The first orbit is dominated by winding of the magnetic field lines and nonlinear growth of the 
MRI. Both of these cause rapid redistributions of disk material and angular momentum. The initial torus 
is stretched radially and material begins to accrete onto the hole and is also carried out to large radii. A 
strong current sheet forms in the initial symmetry plane of the disk through differential winding. 

From orbits 1-2, MRI driven turbulence begins to grow in the inner parts of the disk. At the same time, 
some bending of the disk due to the differential precession from the hole becomes apparent. The MRI is 
fully developed through most of the disk around orbit 2. 
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Fig. 2. — Volume visualization of the logarithm of density (scaled from 0.008p max ,o to 0.8/5 max ,o) a ^ ( a ) * = 0, 
(6) 1, (c) 2, (d) 4, (e) 7, and (/) 10i or b- Half of the disk has been cut away to reveal the cross section. The 
black hole spin axis is oriented vertically in each frame so that the initial torus is tilted 15° to the right. 
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By about orbit 7-8, the disk has reached a quasi-steady state. In the remainder of this section we detail 
the properties of the resultant structure. We follow an "inside-out" track, starting from key features of the 
flow near the hole and working toward larger radii. Where practical, we draw attention to similarities and 
di fferences between the qu asi-steady structure that results in this simulation and the untilted simulations 
of iDe Villiers et al.l (|2003l ). In particular, we draw attention to the fact that some features, such as the 
inner torus and plunging region, are significantly altered, while others, such as the main body and coronal 
envelope, show very similar properties. Again, because of the varying levels of refinement along the poles, 
we do not discuss the evacuated funnel or funnel- wall jet in this paper. 



3.1. Global Structure 

3.1.1. Plunging Streams 

Perhaps the most striking feature in the tilted disk at late times are the two opposing strea ms that start 



from high latitudes both with respect to the black-hole symmetry plane and the disk midplane (jFragile et al 



20071 ) . Figure [3] shows a zoomed-in view of the region around the black hole including these streams. Note 
that stream 1 remains entirely above the black-hole symmetry plane, while stream 2 remains below. Clearly 
the material in each stream is in a plunging orbit into the black hole. Hence, we refer to these features as 
the "plunging streams." 

Figure [4] captures the plunging streams from a different perspective. This image is a view looking down 
the angular momentum axis of the black hole onto a single isodensity surface. The two opposing streams 
are clearly visible in the interior region of the disk as well as two relatively evacuated lobes. 

As material passes through the plunging streams it undergoes strong differential precession. As we show 
below, the precession totals approximately 180°, accounting for how the material in the plunging streams is 
able to enter the black hole from the opposite azimuth from which it began its plunge without ever passing 
through the symmetry plane of the hole. 

Two very important points to make about these streams is that they appear to be stable and stationary. 
They begin forming as early as t = 7t OT h and last until the end of the simulation. During this time their 
azimuthal location does not change appreciably. The interesting questions are why do these opposing plunging 
streams form and why do they start from such high latitude with respect to the black-hole symmetry plane 
and disk midplane? The answers, of course, are related and the fundamental cause is the aspherical nature 
of the gravitational spacetime around the rotating black hole. This is best illustrated by considering the 
dependence of nsco on inclination for orbits that are circular in the sense that they have constant coordinate 
radius. Briefly, nsco is the radius at which the quantity 

dr x 2 



R ee A 2 i^—j = [E{r 2 +a 2 )- aL z ] - A [r 2 + (L z - aE) 2 + Q] (22) 

and its first two derivatives equal zero, i.e. R = R' = R" = 0, where E, L z , and Q are t he energy, an gular 
momentum, and Carter constant, respectively, describin g orbits around Kerr black holes ( Hughes! 200 1[ ) and 



A = r 2 + a 2 cos 2 9 and A = r 2 — 2Mr + a 2 . Following iHughesI (|200lh . we can eliminate Q in favor of the 
inclination i defined as 

C0SZ= (I^%^' (23) 
Figure [5] illustrates this dependence for a few selected cases of a. The key point of the formula and the plot 
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Fig. 3. — Zoomed in view of the inner IOtg of the accretion flow revealing two opposing, high-latitude 
streams of material connecting the disk to the horizon (indicated by arrows). Data is taken from the last 
frame of the simulation (t = 10t ol -b). To emphasize the plunging streams, the scaling in this figure is adjusted 
from that used in Fig. H]by adding a density isosurface at p = 0.024 j o max .o. The figure is oriented as in Fig. [5] 
with the black- hole spin axis vertical. The black-hole symmetry plane (black line) and initial disk midplanc 
(blue line) are marked for reference. Note that stream 1 remains entirely above both planes while stream 2 
remains below. 



-12- 




Fig. 4. — Isodensity contour at p = 0.1p ma x,o from the same time slice as Fig. [3] (t = 10t or b) viewed down 
the angular momentum axis of the black hole. The initial disk angular momentum axis (and polar axis of the 
grid) is tilted 15° to the right in this image. One plunging stream (indicated by solid arrow) starts near the 
left edge of the figure and connects to the hole on the right. This stream lies entirely above the black-hole 
symmetry plane and corresponds to stream 1 in Fig. [3l The opposing stream (stream 2) remains below the 
black- hole symmetry plane and is seen connecting with the horizon on the left. 
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is that orbital stability around a rotating black hole is strongly dependent on the inclination of the orbit. 
Notice that the unstable region increases monotonically for increasing inclination. 

We can make better use of the information in Figure [5] by converting it to a polar plot (using only the 
prograde orbits) and overlaying it onto a plot of data from the simulation, as is done in Figure [6] Such a 
polar plot creates a representation of the prograde "ISCO surface" (symmetric about the spin axis of the 
black hole), which gives a clear indication of where the most unstable regions of the spacetime are. Note 
that the plunging orbits highlighted previously start near where the disk first encounters the ISCO surface. 
More precisely, the streams start near the largest cylindrical radius (rcos$) of the ISCO surface, measured 
with respect to the angular momentum axis of the disk. This explains why the plunging streams start at 
such high inclinations relative to the black-hole symmetry plane and the disk midplane and why there are 
only two streams. The plunging region is no longer azimuthally symmetric from the perspective of the disk. 

Another point to take away from Figures [5] and [5] is that nsco is larger for larger inclinations. Thus, 
for a given black-hole spin, plunging orbits will always start further away from the hole for more tilted disks. 
The tilted black hole effectively acts like an untilted black hole of lower spin, which would likewise have a 
larger nsco 



3.1.2. Inner Torus 

In our tilted simulation, the plunging streams appear to connect directly to the main disk body without a 
clearly identifiable intermediate "inner torus" . This appears to be a particular result of the tilted simulation 
and not, for i nstance, due to the diffe rences in the coordinates used in our simulation (Kerr-Schild) versus 
those used in be Villiers etaD (|2003h (Boyer-Lindquist) or numerical techniques. We base this statement 



on the fact that our ow n untilted simula t ion in Kerr-Schild coordinates shows an inner torus very similar 
to the one described in be Villiers et alj |2003h . For instance, Figure [7| shows the shell-averaged density 



and pressure as a function of radius for our tilted and untilted simulations. Shell averaged quantities are 
computed over the most refined grid as follows: 



o 



(Q) A (r,t) = j I I QV^gdMip, (24) 



where A = ^/—gd-ddcp is the surface area of the shell. The data in Figure [7] has also been time- 

averaged over the final orbit, 9t or b = £ m in < t < i max = 10t or b, where time averages are defined as 



(Q)t = 1 —— / Qdt . (25) 

'■max ^min J t m i a 

In the untilted simulation, both the density and the pressure show local maxima near 4.5rc, indicating an 
inner torus. The tilted simulation, on the other hand, shows only marginal evidence for local maxima near 
10r G . 

Another check of the presence of an inner torus is to look at the distribution of specific angular mo- 
mentum in the disk. Because the inner torus is partially supported by pressure gradients, some portion 
of the flow must be locally super-geodesic. In Figure [8] we plot the density- weighted shell average of the 
specific angular momentum (£)a — (p^}a/(p)a as a function of radius, again time-averaged over the interval 
t = 9 to 10i or b- We compare this against the specific angular momentum distribution of circular orbits with 
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Fig. 5.— Plot of the inclination dependence of rjsco f° r black-hole spins a = 0, 0.5, 0.9, and 0.998. 
Inclinations < i < 90° represent prograde orbits, whereas inclinations 90° < i < 180° represent retrograde 
orbits. 
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Z-Axiz 




X-Axi z 



Fig. 6. — Meridional plot (ip = 0) through the final dump (t — 10i O rb) of the simulation showing a pseudocolor 
representation of the logarithm of density (scaled from 0.008p max .o to 0.8p maXi o as in previous figures) and 
an isocontour of density at p = 0.024p maX: o (red curve). Unlike previous figures, this one is shown oriented 
in the sense of the grid, so that the black hole is tilted 15° to the left. The plot is overlaid with a polar 
plot of the "ISCO surface" for prograde orbits about an a = 0.9 black hole (white curve). This surface 
is symmetric about the spin axis of the hole. Notice that the plunging streams from Figs. [3] and [4] start 
near the largest cylindrical radius (r cos i?) of this surface (indicated by white arrows) and connect with the 
horizon approximately 180° away in azimuth (indicated by black arrows). 
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inclinations of 15° and 0°. These are calculated from the following expression 

Nx + A(Mr) 1 / 2 iV 2 1/2 cosi 



D 



(26) 



where 



Ni = -aMr (3r 2 + a 2 - 4Mr) cos 2 i , (27) 

N 2 = r 4 + a 2 sin 2 i (a 2 + 2r 2 - 4Mr) , (28) 

and 

D = a 2 (2r 2 + a 2 - 3Mr) sin 2 i + r 4 + 4M 2 r 2 - 4r 3 M - Mrs 2 , (29) 

which comes from noting that for circular orbits R = R' = from equation (|22p and from the definition 
£ = L z /E. Both simulations show a nearly geodesic angular momentum distribution through most of the 
disk with a small region of super-geodesic flow inside 10rc- This region clearly corresponds to the inner 
torus in the untilted simulation. It also suggests that there should be an inner torus in the tilted simulation, 
though, again, this is not as evident in the plots of density and pressure. 

Another indication that the inner torus is less prominent in the tilted simulation than the untilted one 
comes from comparing the total rest mass in the near-hole region (r < r cut — IOtq). This is done in Figure 
O where we plot the time histories of the total (volume-integrated) rest mass 

/•2tt /-it rr cut 

(pu°) v = / / DdrdMtp . (30) 

JO JO Jr min 

At t = 10i or b, the inner torus is 42% less massive in Model 915h. 

When present, the inner torus usually performs two functions: regulating the accretion of matter into 
the black hole and serving as the launching point for the funnel-wall jet. Therefore, we may expect a weaker 
funnel-wall jet (to be discussed in future work) and a higher mass accretion rate in our tilted-disk simulation 
relative to the untilted simulation due to the less prominent inner torus in the former. We compute the mass 
accretion rate 

/•27T />7T 

M(r) = / DV r dMcp (31) 
Jo Jo 

100 times per i or b (about every 8M) at each of the external grid boundaries. Figure [TUb shows a plot 
comparing M (r m i n ) for our equivalent tilted and untilted simulations. When averaged over the quasi-steady 
state of each simulation (from t = 7 to 10£ or b), (M) t into the hole for the tilted simulation (915h) is 
7.2 x 10~ 6 , while for the untilted one (90h), it is 4.9 x 10 6 . There is a clear tendancy toward a higher M 
in the tilted-disk simulation. 

Figure \Wb compares M of the tilted disk simulation at three different resolutions. Due to the chaotic 
nature of the mass accretion we do not expect the individual peaks to match; yet we are encouraged that 
the overall shape and magnitude of the two high-resolution models (915h and 915vh) are very consistent, 
suggesting we are reasonably well converged. The medium resolution simulation (Model 915m), on the other 
hand, is clearly underresolved. 



3.1.3. Main Disk Body & Coronal Envelope 



The main disk body does not differ substantially between the tilted and untilted simulations, except in 
the notable fact that the tilted disk precesses (as discussed in §3.2.21 below). Likewise, the coronal envelope, 
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Fig. 7. — Plot of {{p)A)t {solid line) and ((P)a)* [dashed line) as a function of radius for equivalent (a) 
tilted Po = 15° (915h) and (6) untilted (3q = 0° (90h) simulations. For both simulations, the data has been 
time-averaged over the interval t = 9 to 10i O rb- The density and pressure have been normalized by their 
respective maxima at t = 0, which are the same in both simulations. 
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Fig. 8. — Plot of the density-weighted time- and shell-averaged specific angular momentum {{t)A)t [thick 
line) as a function of radius for equivalent (a) tilted (3q = 15° (915h) and (6) untilted /3 = 0° (90h) 
simulations. For both simulations, the data has been time-averaged over the interval t = 9 to 10i O rb- In each 
plot a comparison is provided with the specific angular momentum of circular orbits with inclinations of 15° 
and 0°, respectively [dashed line). For reference we also include the initial angular momentum distribution 
in the midplane of the torus [thin line). 




Fig. 9. — Total rest mass in the near-hole region (r < lOrc) as a function of time for the tilted (915h & 
915vh) and untilted (90h) simulations. The mass and time are normalized by the initial mass and orbital 
period of the torus, respectively. 
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which extends above and below the disk, shows very similar properties in all our simulations. The material 
in the coronal envelope is characterized by low density and rough magnetic equipartition (/3 mag ~ 1)- By 
contrast the main body of the disk is generally gas-pressure dominated (8 mas <C 1). Therefore, a plot 
of 8 mag and p, such as F igure 111! provides a convenient means to identify these two regions. As found 
in |Pe Villiers et al.l (|2003l ). the material in the coronal envelope moves mostly radially outward, yet has 
(—hu t < 1). This suggests that the material may be gravitationally bound, in which case it must circulate 
back to the disk at large radii. However, we point out that this definition of binding energy ignores the 
contribution of the magnetic fields, so some of this material may in fact escape the system. We plan to 
examine outflows from tilted disks more thoroughly in future work. 

Because the disk is precessing, its angular momentum axis does not remain aligned with the grid. 
Therefore, an azimuthal slice through the disk at late times, such as Figure [TTJ may give the impression that 
the disk has aligned with the symmetry plane of the black hole when indeed this is not the case. We now 
turn to the question of disk alignment and precession. 



3.2. Results Specific to A Tilted Disk 

3.2.1. Tilt 

One key diagnostic for describing the global response of a tilted disk subject to Lense-Thirring precession 
is the tilt between the angular momenta of the black hole and disk as a function of radius and time. For 
example, in the Bardeen-Petterson solution, no time variability is observed, and the tilt transitions from 
nearly zero close to the black hole to a non-zero asymptote at large radii. 

As in Fr agile fc Anninosl |2005h . we recover the tilt from the simulation data using the definition 



8{r) 

where 



Jbh • JDisk(r) 



|JBH||JDisk(r)| 



(32) 



Jbh = (-oM sin 8qx, 0, aM cos 8qz) (33) 
is the angular momentum vector of the black hole and 

JDiskM = [(JDisk)li, (JDisk)22/, (</Disk)3-S] (34) 

is the angular momentum vector of the disk in an asymptotically flat space. This is given by 

( Joisk)p = v-g°g a ' (35) 

where 

If = J (i'T" — x u T^°) d 3 x, (36) 

and S a — J T a0 d 3 x. The equations for and S a are integrated over concentric radial shells of the 
most-refined grid layer, e.g. 

crO 



S a (r) = / / T^V^Ardtfdp . (37) 

The unit vector y points along the axis about which the black hole is initially tilted and z points along the 
initial angular momentum axis of the disk. 
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In Figure [12l we show the radial profile of (3 time averaged over the interval 9t rb <t< 10i O rb- Recall 
(3q = 15° for this simulation. This profile remains fairly consistent over many orbital times once the quasi- 
steady state is reached, so the time-averaged data gives a good representation for all t > 7t rb- The variability 
from this time-averaged profile is generally < 20% and is generally carried by moderate amplitude waves 
traveling through the disk. The increase in tilt at r < lOr^ is attributable to the high latitude plunging 
streams described in ^3. 1.11 

One very obvious characteristic of the profile in Figurc[T2"lis that (3 does not approach zero except perhaps 
very close to the hole. Thus we do not see evidence for the Bardeen-Petterson effect in this simulation. This 
is not surprising since the Bardeen-Petterson solution is only expected for thin disks (H/r < a). This is not 
the applicable regime for this simulation, as we illustrate in Figure Q31 which shows H/r and a plotted as 
functions of r. The scale height H(r) is defined in each radial shell as one-half the distance (0.5r A$) between 
the two points where p = pmax/e, where we use the time-averaged density along the tp = azimuthal slice. 
The dimensionless stress parameter a in the disk is taken to be 



a 



\u r uV\\B\\ 2 - B r B<?\ 
4nP 



We restrict the calculation of a to only bound material (—hut < 1). Using these definitions we find H/r 
and a < 0.01 through most of the disk. 



(38) 



0.2 



Since warps in slim disks are expected to propagate as ben ding waves, it may se em unusual at first that 
we see little evidence for such waves in Figure [T^J For instance, iLubow et al.l (|2002l ) provides an analysis of 
the theory of bending waves in nearly Keplerian, weakly inclin ed disks and predicts tha t the tilt (3 should 
be a time-indep endent, oscillatory f unction of radius (see also Markovic k, Lambl 1998 ). However, using 
equation (16) of ILubow et al.l (J2j302j), we estimate the wavelength of such oscillations for our simulation to 
be 

^ -■ (39) 



— - 50M 



A ~ {6a) 1 / 2 

at r = lOrg. This is strongly radially dependent (A oc r 9 / 4 with H/r ~ constant), so oscillations of (3 are 
essentially absent outside r = lOr^, consistent with what is shown in Figure [T2l 

The same conclusion, that f3 is no t expected to oscilla te outside r = lOrp for this simulation, is also 
reached by considering equation (22) of Lubow et al. ( 2002f h That equation defines a dimensionless variable 



24a 

,2 



1/2 r -(h+l/4) 

h + 1/4 



(40) 



which is used to identify the transition radius between oscillatory behavior and asymptotic behavior, where 
h and e are used to parameterize the radial dependence of the disk scale height H/r — er' 1-1 . Whenever 
x » 1 (small r), oscillations should be prominent, whereas whenever x « 1 (large r), (3 tends to the 
outer boundary value. For our simulation, with e » 0.2 and fowl, 1 = 1 at r « 10rc- Thus, from both 
approaches, it is clear that our simulation does not satisfy the criteria to develop large oscillations in (3 
within the main body of the disk. 

Inside r = lOrc the density of the disk drops off rapi dly and the dynam ics are dominated by the 
plunging streams, which are not accounted for in the model of Lubow et al. ( 2002h . Nevertheless, we appear 
to capture one-half of one wavelength of a bending wave oscillation inside r = lOry;, based on Fig ure [T2l 
Thus, overall our results seem to be generally consistent with the predictions of ILubow et al.l ( 20021) . 
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3.2.2. Precession 

A second useful diagnostic for tilted disks is the twist 7 of the disk as a function of radius and time. 
We define the precession angle (twist) as 



/ \ [ J BH X JDisk(f) .1 / .-, s 

7(r) = arccos — — — • y , (41) 

_ |JbH X JDi s k(r)| 

From this definition, 7(7-) = throughout the disk at t = 0. In order to capture twists larger than 180°, 
we also track the projection of Jbh x JDisk(^) onto x, allowing us to break the degeneracy in arccos. A 
time-averaged plot of 7 is provided in Figure 1141 

As described in our previous work ([Fragile fc Annino sl l2005h . we expect differential Lense-Thirring pre- 



cession to dominate w henever the precession timescale tyr = ^it = Q u /q is shorter than local dynamical 



timescales in the disk (jBardeen fc Pettersonlll975tlKumar fc Pringldll985f ). We consider three possible limit- 
ing timescales: the mass accretion timescale i a cc = r/V T , where V = ((pV t )a/ {p)A)t is the density-weighted 
average inflow velocity; the sound-crossing time t cs = r/c s , where c s = ((pc s )a/ (p)A)t is a density-weighted 
average of the local sound speed; and the Alfven crossing time Ia — r/V a, where V a is a density- weighted 
average of the local Alfven speed. The local sound speed is recovered from the fluid state through the relation 



= T(r - l)P/[(r - l)p + TP}. The Alfven speed is 



\B\ 



VA 



4TT P h+ \\B\\ 2 



(42) 



Since c s and va are defined in the frame of the fluid, it is not strictly accurate to compare t cs and tA to 
quantities defined using the coordinate time (such as £lt and f2 _1 ). However, we are mostly concerned with 
the timescales in the main body of the disk where such discrepancies are small. From Figure HSl we can see 
that the Lense-Thirring precession timescale is longer than the sound-crossing time at virtually all radii. 

Since the sound-crossing time is short compared to the precession timescale throughout the bulk of the 
disk, pressure waves strongly couple the disk material. The disk, thus, responds as a single entity to the 
torque of the black-hole and precesses as a global structure. Such global precession has been n oted before in 
low Mach number hydrodynamic disks ( Nelson fc Papaloizou 20001 Fragile fc Anninos 20051 ). To estimate 
the precession period, we have plotted 7, averaged over the bulk of the disk (20 < r/rc < 50), as a function 
of time in Figure [TBI A linear fit to this plot yields a precession period of T prcc w 0.3(M/M Q ) s, which 
corresponds to about 80t or b- This is longer than the evolution time of all of our models, so we have had 
to extrapolate the full precession period. However, Model 915m is run to 20t O rb and shows a nearly linear 
growth of precession over the full simulation. 

Classically, we expect the precession perio d for a solid-body r otator with angular momentum J subject 
to a torque r to be T proc = 27r(sin/3) (J/r) (|Liu fc Melial l2002n . Assuming a radial dependence to the 
surface density of the form E = Ei(r/ri) _l » and ignoring higher order general relativistic corrections, we have 
J = 27rAf 1 /2 Sir C r .5/2-C [1 _ (r . i/ro) 5/2-C]/(5/ 2 _ C ) and T = ATT(sinp)aM 3 / 2 ^[l-(r t /r o ) 1 / 2 +<}/[r 1 l /2 (l/2+0}, 
where and r Q are the inner and outer radii of the evolved disk, respectively. Therefore, 

(r l /r ) 5 / 2 -q 



T 



prcc 



^ (1+2C)r ^ r l/2 + C {1 

(5 - 2C) aM [1 - {n/ro) 1 '^} 



(43) 



For Ti — 10r<3, r Q — 50rc, and £ = (the value we find in our simulation), equation (|43[) predicts T prcc = 
0.3(M/M Q ) s, which is the same as the observed value in the simulation. Note that equation (|4"3"| differs 
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from the test particle Lense-Thirring precession period because T prec depends on the total torque integrated 
over the entire disk. 



Discussion 



In this paper we studied the evolution of an MRI turbulent disk that was tilted with respect to the 
spin axis of a modestly fast rotating black hole. Although this prescription can lead to a Bardeen-Petterson 
configuration for some disk parameters, we did not see evidence for this in this simulation, as alignment of 
the disk with the equatorial plane of the black hole did not occur. This is not surprising since this simulation 
was ca rried out in the thick-d isk regime where H/r > a and warps produced in the disk propagate as 
waves (|Papaloizou fc Lin 1995 ) , rather t han diffusively as in the Bardeen-Petterson case. Since the expected 
bending wavelength |Lubow et al" 2002h turned out to be longer than the radial extent of the disk in the 
simulation, little warping of the disk was observed. Instead the unwarped disk precessed uniformly. The 
extrapolated precession period T plec w 0.3(M/M Q ) s equates to periods of « 3 s and « 3 d for black holes of 
mass M — 10M & and M — 10 6 M Q , respectively. Such global disk precession could explain certain variability 
features observed from accreting black ho les, such as low-frequency QPOs (LFQPOs) (jStella et al.lll999t 
Liu fc Melia 20021 : Schnittman et al. 2006 ). since the observer's viewing angle of the inner, X-ray emitting 
region of the disk would vary periodically. 

If the inner disk is optically thick enough to produce relativistically-broadened reflection features, such 
as an iron Ka line, then such precession should a lso be observable through periodic changes in both the 
shape and strength of the lines (jFragile et al.l 120051 ). These changes should be corr elated with the phase o f 
the corresponding LFQPO. Such a correlation has been observed in GRS 1915+105 (jMiller fc Hom an 2005), 
although only between line strength and QPO phase; those data were not sufficiently resolved to determine 
the line shape. 



Generally, we expect the precession period to be given by equation (j43|) , which has a strong dependence 
on the radial distribution of the disk (oc r'J 2 ^r*^ 2+ ^). One idea to consider is that the outer radius 
may correspond to the truncation r adius proposed to e xplain the hard state of black hole X-ray binaries 



(e.g. 



Esin et al.l (| 19971) . but see also iRvkoff et al.l (|2007l) ). In this case our simulated disk would represent 
the hot, geometrically thick flow that fills the region inside the truncation radius. The LFQPO would 
then correspond to t he precession frequency of this inner flow, in which case it should scale as r 5 / 2+ ^ . 
Sobczak et alj (|2000h explored the dependence of the LFQPO frequency on spectral fitting parameters, 
including what would be the truncation radius in the context of the suggested hard state model. They 
studied two sources, XTE J1550-564 and GRO J1655-40, and found opposite trends between frequency and 
radius. For XTE J1550-564 the observed frequency was ^lfqpo ~ 5 Hz, and the observed truncation radius 
was r /r G = 2.7(WM Q /M)(D/6 kpc)(cos #) -1 / 2 . From equation <[43j) we would expect 



2C 



tt(1+2C) 



2/(5-0 



(-) 



2/(5-0 



I'G 



-(1+20/(5-20 



(uM) 



-2/(5-0 



(44) 



In our simulation we found C ~ 0, which gives r a « 33r<3 for M = IOMq and v = 5 Hz. This is considerably 
larger than the observed value. However, some of the discrepancy may be attributable to the large uncer- 
tainties in the parameters used to describe this source, including its distance, mass, and inclination. Also, if 
the surface density in XTE J 1550-564 depends strongly on radius, which was not the case for our simulated 
disk, then our prediction would change significantly. Further observational studies along this line are needed 
to test this prediction more thoroughly. 
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Although the main body of the disk was not significantly altered by the tilt, we did find significant 
differences in the inner regions of the flow when compared with untilted simulations. First, a tilted disk 
encounters the generalized ISCO surface at a larger radius than an untilted disk. This causes the plunging 
region to start further out. The binding energy of the innermost material in the disk is therefore less than 
it would be for an aligned disk, and the overall radiative efficiency should then be reduced. 

On the other hand, tilting the disk appears to produce a higher overall mass accretion rate (shown 



here in Figure 10a; also discussed in lLodato &: Pringld I2006T ) . A tilted accretion disk will therefore have 



a lower surface density than an untilted disk with the same accretion rate. This may affect the emergent 
spectr um, especial l y for hot, optically thin flows. On the other hand for flows that are effectively optically 



thick, IDavis et alJ (|2005l ) found that the emergent spectra are remarkably independent of the overall stress 



and surface density. 

We also found that the plunging region is not axially symmetric. Instead, accretion onto the hole in the 
tilted-disk case occurs through two discrete streams of material that leave the disk at high latitudes with 
respect to the black-hole and disk symmetry planes. This may affect the magnitude of magnetic torques 
exerted by the plunging region on the disk. An interesting question for future work is how these streams vary 
on the timescale of the precession of the disk. We intend to explore the detailed properties of the plunging 
region and innermost disk in a future paper. 

The tilted disk also seems not to have formed a clearly identifiable inner torus. This could be significant 
because the inner torus serves as a launching point for the matter-dominated, funnel-wall jet. The absence 
of a prominent inner torus may lead to a weaker matter jet. However, the present simulation is not suited 
to addressing this issue because of the poor and varying resolution used near the pole. Instead, we plan to 
explore jets and outflows from tilted disks in future work. 

In many respects the tilted disk simulation exhibited properties consistent with an untilted disk around 
a black hole of lower spin. These included the larger plunging radius, higher mass accretion rate, and 
less prominent inner torus. Thus black-hole tilt could hamper efforts to estimate black-hole spin based on 
such properties. Indeed, it is commonly stated that astrophysical black hole spacetimes depend on just 
two parameters: mass and spin. But it should be remembered that the observed properties of black hole 
accretion disks also depend on their inclinations with respect to the spin axes of their central black holes. 
This inclination should be a target of future observational programs that use accretion disks as surrogates 
to study properties of black holes. 
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Fig. 10. — (a) Plot of the mass accretion history from Model 915h with flo — 15° (thick line) and Model 90h 
with (3 = 0° (thin line). The accretion rate and time are normalized by the initial mass and orbital period 
of the torus, respectively, (b) Plot of mass accretion rate, comparing our medium (915m), high (915h), and 
very high (915vh) resolution tilted disk simulations. The very high resolution simulation was only run to 
t = 5i or b- 
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Fig. 11. — Azimuthal slice through the simulation along ip = taken from the final dump (t = 10i O rb)- 
The ratio of magnetic pressure to gas pressure (/3^a g ) is represented as a pseudocolor plot. The colors are 
scaled logarithmically and cover the range 10~ 2 < /3ma g < 10 2 . The gas density is given by isocontours at 
p = 10~ 2 , 10~ 1,5 , 10 , and 10~ a5 p max ,o. As with Fig. [HI this figure is oriented in the sense of the grid, 
so that the black hole is tilted 15° to the left. The apparent tilt of the disk is actually due its precession 
about the black-hole spin axis, such that the angular momentum axis of the disk is no longer in the plane 
of this image; the disk has not actually realigned with the hole. We remind the reader that the region near 
the poles is not sufficiently resolved, so caution should be used when interpreting results there. 
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Fig. 12. — Plot of the tilt {/3)t as a function of radius through the disk. The data for this plot has been time 
averaged from t = 9 to 10i O rb- The initial tilt was @o = 15°. 
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Fig. 13. — Plot of the scale height (H) t /r and magnetic stress parameter (a) t , time averaged over the interval 
7i rb < t < lOiorb- This plot shows that this simulation falls into the thick-disk limit H/r > a. 
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Fig. 14. — Plot of the twist (7)4 as a function of radius through the disk. The data for this plot has been 
time averaged from t — 9 to 10t or b- Initially the twist was zero throughout the disk. The disk matter has 
precessed roughly <~ 180° by the time it reaches the hole. The shape of this twist profile remains fairly 
constant throughout the simulation. 
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Fig. 15. — Plot comparing various timescales within the disk, including the Lcnse-Thirring precession 
timescale iur, the accretion timescale i acc , the sound-crossing time t cs , and the Alfven crossing time t^. 
All timescales are normalized by the local orbital period in the midplane of the black hole, The data 

for this plot has been time averaged from t — 9t or b to t — 10i or b- 
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Fig. 16. — Plot of the twist 7, averaged over the bulk of the disk (20 < r/rc < 50), as a function of time. The 
slope of this plot can be used to estimate the precession period of the disk as a whole, which is 0.3(M/M Q ) 
s. 



